Source of this introduction: Cartopy Documentation
First, we have to import Matplotlib and Cartopy including its Coordinate Reference System.
## Import modules
# Matplotlib
%matplotlib notebook
import matplotlib.pyplot as plt
plt.ioff() # turn off interactive plotting mode
# Cartopy
import cartopy as cp
# Coordinate Reference System
import cartopy.crs as ccrs
Bad key "axes.color_cycle" on line 17 in /Users/awi/.matplotlib/stylelib/mystyle.mplstyle. You probably need to get an updated matplotlibrc file from http://github.com/matplotlib/matplotlib/blob/master/matplotlibrc.template or from the matplotlib source distribution
We can access the individual projections from the imported cartopy coordinate reference system, e.g. ccrs.PlateCarree().
plt.figure(figsize=(10,5))
# Instantiate matplotlib axes with Cartopy Plate Carrée projection
ax = plt.axes(projection=ccrs.PlateCarree())
# Add coastlines
ax.coastlines();
plt.show()
Displaying available projections
Explore Cartopy's projections list or check
dir(ccrs)
# or hit TAB in an IPython Console
ccrs.<TAB>
directly in your console!
# Check for available projections
print dir(ccrs)
['ABCMeta', 'AlbersEqualArea', 'AzimuthalEquidistant', 'CRS', 'EuroPP', 'GOOGLE_MERCATOR', 'Geocentric', 'Geodetic', 'Geostationary', 'Globe', 'Gnomonic', 'InterruptedGoodeHomolosine', 'LambertAzimuthalEqualArea', 'LambertConformal', 'LambertCylindrical', 'Mercator', 'Miller', 'Mollweide', 'NearsidePerspective', 'NorthPolarStereo', 'OSGB', 'OSNI', 'Orthographic', 'PROJ4_VERSION', 'PlateCarree', 'Projection', 'Robinson', 'RotatedGeodetic', 'RotatedPole', 'Sinusoidal', 'SouthPolarStereo', 'Stereographic', 'TransverseMercator', 'UTM', 'WGS84_SEMIMAJOR_AXIS', 'WGS84_SEMIMINOR_AXIS', '_BoundaryPoint', '_CylindricalProjection', '_RectangularProjection', '_Satellite', '_WarpedRectangularProjection', '__builtins__', '__doc__', '__document_these__', '__file__', '__name__', '__package__', '__warningregistry__', '_ellipse_boundary', '_find_first_ge', 'absolute_import', 'abstractproperty', 'cartopy', 'division', 'epsg', 'math', 'np', 'prep', 'print_function', 'sgeom', 'six', 'warnings']
Let's choose an equal-area projection, e.g. Mollweide
# Mollweide projection with central longitude 180° W
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=180))
ax.coastlines()
# Add stock background image
ax.stock_img()
plt.show()
Let's have a look at another one: Interrupted Goode homolosine projection
# Intantiate new figure
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.InterruptedGoodeHomolosine())
ax.stock_img()
plt.show()
Cartopy feature interface relys on the comprehensive collection of free vector and raster map data provided by Natural Earth Data.
Natural Earth Data provides at 1:10m, 1:50m, and 1:110m scales.
Most commom features are pre-defined:
| Feature | Â Description |
|---|---|
| cartopy.feature.BORDERS | Country boundaries |
| cartopy.feature.COASTLINE | Coastline, including major islands |
| cartopy.feature.LAKES | Natural and artificial lakes |
| cartopy.feature.LAND | Land polygons, including major islands |
| cartopy.feature.OCEAN | Ocean polygons |
| cartopy.feature.RIVERS | Single-line drainages, including lake centerlines |
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
# Add oceans and high resolution shapefile of lakes
ax.add_feature(cp.feature.OCEAN)
ax.add_feature(cp.feature.LAKES.with_scale('10m'))
plt.show()
Let's integrate our custom feature.
Example: I want to plot railroads.
Get download link of the data provided by Natural Earth Data: http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_railroads.zip
Figure out category (='cultural'), name (='railroads'), and scale (='10m').
Define your railroads feature (Cartopy gets the data for you):
railroads = cp.feature.NaturalEarthFeature(category='cultural',
name='railroads',
scale='10m',
facecolor='none', # facecolor is set to 'none', because the polygons should not be color-filled
color='blue') # color of the railroads
ax.add_feature(railroads)# Railroads feature example
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
## Add a custom feature from Natural Earth Data: world's railroads
railroads = cp.feature.NaturalEarthFeature(category='cultural',
name='railroads',
scale='10m',
facecolor='none', # facecolor is set to 'none', because the polygons should not bet filled
color='blue')
ax.add_feature(railroads)
<cartopy.mpl.feature_artist.FeatureArtist at 0x1809f56fd0>
plt.show()
Data are added to a projection just the same as they are added to a standard matplotlib figure.
Example: Plotting airlines between several cities
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
# Hamburg coordinates
hh_lon, hh_lat = 10, 53
# Portland coordinates
pl_lon, pl_lat = -122, 45
# Sydney coordinates
sd_lon, sd_lat = 151, -33
# Plot connecting lines on a plane
plt.plot([pl_lon, hh_lon], [pl_lat, hh_lat],
color='gray', linestyle='--', linewidth=3,)
plt.plot([sd_lon, hh_lon], [sd_lat, hh_lat],
color='gray', linestyle='--', linewidth=3)
# Add labels for each cities
plt.text(pl_lon - 3, pl_lat - 8, 'Portland',
horizontalalignment='right', fontsize=18)
plt.text(sd_lon - 3, sd_lat - 8, 'Sydney',
horizontalalignment='right', fontsize=18)
plt.text(hh_lon + 3, hh_lat + 8, 'Hamburg',
horizontalalignment='left', fontsize=18)
Text(13,61,'Hamburg')
plt.show()
Geocentric
Cartersian coordinate system, where x, y, z are coordinates from the center of the Earth
Geodetic
Latitude/longitude coordinate system with spherical topology, geographical distance and coordinates in degrees
To get the actual airlines between two cities on Earth, the data has be transformed to spherical coordinates using the keyword transform and the method ccrs.Geodetic:
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
# Plot connecting lines on a plane
plt.plot([pl_lon, hh_lon], [pl_lat, hh_lat],
color='gray', linestyle='--', linewidth=3,)
plt.plot([sd_lon, hh_lon], [sd_lat, hh_lat],
color='gray', linestyle='--', linewidth=3)
# Add labels for each cities
plt.text(pl_lon - 3, pl_lat - 8, 'Portland',
horizontalalignment='right', fontsize=18)
plt.text(sd_lon - 3, sd_lat - 8, 'Sydney',
horizontalalignment='right', fontsize=18)
plt.text(hh_lon + 3, hh_lat + 8, 'Hamburg',
horizontalalignment='left', fontsize=18)
# Plot connecting lines on spheric coordinates
plt.plot([pl_lon, hh_lon], [pl_lat, hh_lat],
color='blue', linewidth=3, marker='o',
transform=ccrs.Geodetic())
plt.plot([sd_lon, hh_lon], [sd_lat, hh_lat],
color='red', linewidth=3, marker='o',
transform=ccrs.Geodetic())
[<matplotlib.lines.Line2D at 0x1809727e90>]
plt.show()
transform and projection¶Cartopy's core concept is that the projection of your axes is independent of the coordinate system your data is defined in.
The projection argument is used when creating plots and determines the projection of the resulting plot (i.e. what the plot looks like).
The transform argument to plotting functions tells Cartopy what coordinate system your data are defined in.
Cartopy always assumes that the your data are defined in the same coordinate system as the projection you use. So, we have to tell Cartopy the data's coordinate system $\Rightarrow$ most of the time lon/lat grids with cartesian coordinate system (ccrs.PlateCarree).
The safest thing to do is always provide the
transformkeyword regardless of theprojectionyou are using, and avoid letting Cartopy make assumptions about your data’s coordinate system. Doing so allows you to choose any map projection for your plot and allow Cartopy to plot your data where it should be.
Read in netCDF data
# Load netCDF4 module
import netCDF4 as nc4
# NCEP reanalysis data: Monthly 2m air temperature
# On ICDC server:
# filename='/data/icdc/reanalyses/ncep_reanalysis1/DATA/2m_airtemp_monthly/air2m.mon.mean.nc'
filename='Data/air2m.mon.mean.nc'
# Open file and create file identifier
f = nc4.Dataset(filename,'r')
# Retrieve variables
lat = f.variables['lat'][:]
lon = f.variables['lon'][:]
air = f.variables['air'][:,:,:]
# Calculate climatologic mean for July
air_jul = air[6::12,:,:].mean(axis=0)
# Kelvin to Celsius
air_jul = air_jul - 273.15
Add the data on a projection
# Instantiate projection
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=0))
# Plot the data in filled contours; Using transform=ccrs.PlateCarree()
# we tell Cartopy the coordinate system the data is defined in
ax.contourf(lon, lat, air_jul[:,:], transform=ccrs.PlateCarree())
# Add coastlines
ax.coastlines()
# Set the extent of the Axes to the limits of the projection
ax.set_global()
plt.show()
# Import add_cyclic_point methd
from cartopy.util import add_cyclic_point
# Close the gap in the data and the longitude array
air_jul, lon = add_cyclic_point(air_jul, coord=lon)
# Instantiate projection
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Mollweide(central_longitude=0))
# Plot the data in filled contours; Using transform=ccrs.PlateCarree()
# we tell Cartopy the coordinate system the data is defined in
ax.contourf(lon, lat, air_jul[:,:], transform=ccrs.PlateCarree())
# Add coastlines
ax.coastlines()
# Set the extent of the Axes to the limits of the projection
ax.set_global()
plt.show()
Example: 2m air temperature; climatologic mean for July; land without Antarctica
# Numpy
import numpy as np
# Import coordinates formating tool
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
# Instantiate projection
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Miller(central_longitude=0))
# Mask ocean
ax.add_feature(cp.feature.OCEAN, zorder=1) # z order must be higher than the data
# Add coastlines
ax.coastlines()
# Add gridlines
ax.gridlines(color='k', linestyle='dashed')
# Cut off Antarctica
ax.set_extent([-180,180,-60,90]) # lon west, lon east, lat south, lat north
# Plot the data in filled contours, 17 levels from -36 to 36 deg C, specify colormap, extend colormap in both ways
p = ax.contourf(lon, lat, air_jul[:,:],
levels=np.linspace(-36, 36, 17),
extend='both',
transform=ccrs.PlateCarree(),
cmap=plt.get_cmap('RdBu_r'))
# Add colorbar
cbar = plt.colorbar(p, orientation='vertical', label='$^\circ$C', extend='both')
# Add longitude and latitude coordinates
ax.set_xticks([-180, -120, -60, 0, 60, 120, 180], crs=ccrs.PlateCarree())
ax.set_yticks([-60, -30, 0, 30, 60, 90], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# Add title
plt.title('2m air temperature on land -- climatologic mean for July since 1948');
plt.show()
# Read ncep monthly wind data from ICDC
# on ICDC server:
# vwind = '/data/icdc/reanalyses/ncep_reanalysis1/DATA/2m_airtemp_monthly/vwnd10m.mon.mean.nc'
# uwind = '/data/icdc/reanalyses/ncep_reanalysis1/DATA/2m_airtemp_monthly/uwnd10m.mon.mean.nc'
fv = nc4.Dataset('Data/vwnd10m.mon.mean.nc','r')
fu = nc4.Dataset('Data/uwnd10m.mon.mean.nc','r')
# Retrieve variables
lat = fv.variables['lat'][:]
lon = fv.variables['lon'][:]
vwnd = fv.variables['vwnd'][:,:,:]
uwnd = fu.variables['uwnd'][:,:,:]
# Calc average winds for January and July
u_jan = uwnd[0::12,:,:].mean(axis=0)
u_jul = uwnd[6::12,:,:].mean(axis=0)
v_jan = vwnd[0::12,:,:].mean(axis=0)
v_jul = vwnd[6::12,:,:].mean(axis=0)
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Robinson())
ax.stock_img()
# Make streamplot
ax.quiver(lon, lat, u_jan, v_jan, transform=ccrs.PlateCarree(), regrid_shape=(len(lon)/4, len(lat)/4), width=0.001)
plt.title('Winds in January', fontsize=20)
plt.show()
plt.figure(figsize=(10,5))
ax = plt.axes(projection=ccrs.Robinson())
ax.stock_img()
# Make streamplot
ax.quiver(lon, lat, u_jul, v_jul, transform=ccrs.PlateCarree(), regrid_shape=(len(lon)/4, len(lat)/4), width=0.001)
plt.title('Winds in July', fontsize=20)
plt.show()
Of course, these plots need some polishing.
import numpy as np
# Read in netcdf data on global tripolar grid: N deposition
ds = nc4.Dataset('TP04_ndepo_CMIP_NCAR_CCMI-1-0_gr_185001-185012-clim.nc','r')
# Retrieve variables
lat = ds.variables['lat'][:,:]
lon = ds.variables['lon'][:,:]
ndepo = ds.variables['ndepo'][:,:,:]
plt.figure(figsize=(10,3))
ax1 = plt.subplot(121, projection=ccrs.Robinson(central_longitude=0))
ax1.coastlines(resolution='110m', linewidth=0.5)
ax1.set_global()
## Plot using matplotlib contourf
p = ax1.contourf(lon, lat, ndepo[0,:,:],
extend='max',
levels=np.linspace(0.25e-11, 2e-11, 5),
cmap = plt.get_cmap('Reds'),
transform=ccrs.PlateCarree())
cbar = plt.colorbar(p, orientation='horizontal')
cbar.set_label('Ndepo')
ax1.set_title('Plot using contourf')
## Plot using matplotlib pcolormesh
ax2 = plt.subplot(122, projection=ccrs.Robinson(central_longitude=0))
ax2.coastlines(resolution='110m', linewidth=0.5)
ax2.set_global()
p = ax2.pcolormesh(lon, lat, ndepo[0,:,:],
vmin=0.25e-11,
vmax=2e-11,
cmap = plt.get_cmap('Reds'),
transform=ccrs.PlateCarree())
cbar = plt.colorbar(p, orientation='horizontal')
cbar.set_label('Ndepo')
ax2.set_title('Plot using pcolormesh')
Text(0.5,1,'Plot using pcolormesh')
plt.show()
This was only a short introduction into Cartopy - there is a lot more to explore!
Start here:
Choose a feature on Natural Earth Data and visualize it in CartoPy (choose a projection you like).
Make a figure of the course of the river Elbe (or any river you like) - from its source to its outlet.